//
//  Move.cpp
//  phi_four
//
//  Created by Snir Gazit on 3/10/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Move.h"
// General Move

double addRatio(int m, int n) {
	return double(1 + 2 * m) / double(2 * n + NF);
}
double removeRatio(int m, int n) {
	return double(2 * n + NF - 2) / double(2 * m - 1);
}

void Move::init(MTGen* eng, Lattice * lat) {
	eng_ = eng;
	lat_ = lat;
	zeroOne_.init(eng, UniRealParams(0, 1));
}

bool Move::accept(double p) {
	if (p > 1)
		return true;
	double draw = zeroOne_.draw();
	if (p > draw)
		return true;
	else
		return false;
}
bool Move::checkMove() {
	int totalBonds;
	return lat_->testConfiguration(totalBonds);
}

// Shift Move

void ShiftMove::init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
		BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseMove_.init(eng, UniIntParams(0, 2 * 2 * DIM - 1));
//	createRemove_.init(eng, BinomialVariableParams(1, 0.5));
	//openCloseRatio_ = qmcParameters.openCloseRatio.getVal();
	setCoupeling(boltzRatio);
}

void ShiftMove::setCoupeling(BoltzWeightType boltzRatio) {
	RemoveRatio_[0] = 0.;
	for (int j = 0; j < MAX_BOND; j++) {
		CreateRatio_[j] = boltzRatio[j];
		if (j != (MAX_BOND - 1))
			RemoveRatio_[j + 1] = 1. / (boltzRatio[j]);
	}
}

bool ShiftMove::makeMove() {
// Get Ira
	Lattice::SiteIter ira = lat_->ira_;
	int sec = lat_->sectorType_;
	// Choose Move
	int move = chooseMove_.draw();
	// Create Or Remove
	CreateRemove createRemove = (move < 2 * DIM) ? CREATE : REMOVE;
	// Direction
	int dir = move % (2 * DIM);
	Lattice::SiteIter nSite = ira->nSiteArray_[dir];
	Site::BondArray& bondArray = ira->bondArray_;
	BondValType &bondVal = (*bondArray[dir])[sec];
	// Create Move
	if (createRemove == CREATE) {
		// Annihilate Worm
		if (nSite == lat_->masha_) {
			double p = 1. / double(bondVal + 1);
			if (this->accept(p)) {
				bondVal++;
				// lat_->setIra(nSite); (needed??)
				lat_->isClosed_ = true;
				lat_->ira_ = nSite;
				return true;
			}
			// Create New Bond
		} else {
			// Get Total Bonds
			BondValType &totalBondsType = nSite->totalBondsArray_[sec];
			BondValType &totalBondsOther = nSite->totalBonds_;
			double p = (CreateRatio_[totalBondsOther]
					* addRatio(totalBondsType, totalBondsOther))
					/ double(bondVal + 1);
			if (this->accept(p)) {
				bondVal++;
				totalBondsType++;
				totalBondsOther++;
				lat_->ira_ = nSite;
				assert(this->checkMove());
				return true;
			}
		}
	} else if (bondVal != 0) {
		// Annihilate Worm
		if (nSite == lat_->masha_) {
			// Get Total Bonds Ira
			int &totalBondsTypeIra = ira->totalBondsArray_[sec];
			int &totalBondsOtherIra = ira->totalBonds_;
			// Get Total Bonds Masha
			int &totalBondsTypeMasha = nSite->totalBondsArray_[sec];
			int &totalBondsOtherMasha = nSite->totalBonds_;
			double removeIra = RemoveRatio_[totalBondsOtherIra]
					* removeRatio(totalBondsTypeIra, totalBondsOtherIra);
			double removeMasha = RemoveRatio_[totalBondsOtherMasha]
					* removeRatio(totalBondsTypeMasha, totalBondsOtherMasha);
			double p =  removeIra * removeMasha * double(bondVal);
			if (this->accept(p)) {
				bondVal--;
				assert(bondVal < MAX_BOND);
				totalBondsTypeIra--;
				totalBondsTypeMasha--;
				totalBondsOtherMasha--;
				totalBondsOtherIra--;
				lat_->isClosed_ = true;
				lat_->ira_ = nSite;
				// lat_->setIra(nSite); (Needed?)
				return true;
			}
		}
		// Remove Bond
		else {
			// Get Total Bonds
			int &totalBondsType = ira->totalBondsArray_[sec];
			int &totalBondsOther = ira->totalBonds_;
			double p = (RemoveRatio_[totalBondsOther]
					* removeRatio(totalBondsType, totalBondsOther))
					* double(bondVal);
			if (this->accept(p)) {
				bondVal--;
				assert(bondVal < MAX_BOND);
				totalBondsType--;
				totalBondsOther--;
				lat_->ira_ = nSite;
				assert(this->checkMove());
				return true;
			}
		}
	}
	return false;
}

// Create Pair Move

void CreatePairMove::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters_, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
//boltzRatio_ = boltzRatio;
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	chooseMove_.init(eng, UniIntParams(0, 2 * 2 * DIM - 1));
	boltzRatio_ = boltzRatio;
}

bool CreatePairMove::makeMove(int sec) {
	// JumpSite
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	// Move
	int move = chooseMove_.draw();
	// Create Or Remove
	CreateRemove createRemove = (move < 2 * DIM) ? CREATE : REMOVE;
	// Direction
	int dir = move % (2 * DIM);
	Lattice::SiteIter nSite = jump->nSiteArray_[dir];
	BondArray &bondArray = jump->bondArray_;
	int &bondVal = (*bondArray[dir])[sec];
	// Create Move
	if (createRemove == CREATE) {
		// Get Total Bonds
		// Get Total Bonds Ira
		int &totalBondsTypeMasha = jump->totalBondsArray_[sec];
		int &totalBondsOtherMasha = jump->totalBonds_;
		// Get Total Bonds Masha
		int & totalBondsTypeIra = nSite->totalBondsArray_[sec];
		int &totalBondsOtherIra = nSite->totalBonds_;
		double addIra = boltzRatio_[totalBondsOtherIra]
				* addRatio(totalBondsTypeIra, totalBondsOtherIra);
		double addMasha = boltzRatio_[totalBondsOtherMasha]
				* addRatio(totalBondsTypeMasha, totalBondsOtherMasha);
		double p = addIra * addMasha / double(bondVal + 1);
		if (this->accept(p)) {
			bondVal++;
			totalBondsOtherIra++;
			totalBondsTypeIra++;
			totalBondsOtherMasha++;
			totalBondsTypeMasha++;
			lat_->ira_ = nSite;
			lat_->masha_ = jump;
			// lat_->setIra(nSite); (needed??)
			lat_->sectorType_ = sec;
			lat_->isClosed_ = false;
			assert(this->checkMove());
			return true;

		}
	} else {
		double p = double(bondVal);
		if (this->accept(p)) {
			bondVal--;
			lat_->ira_ = nSite;
			lat_->masha_ = jump;
			lat_->sectorType_ = sec;
			lat_->isClosed_ = false;
			assert(this->checkMove());
			return true;

		}
	}
	return false;
}
void CreatePairMove::setCoupeling(BoltzWeightType boltzRatio) {

	boltzRatio_ = boltzRatio;
}

